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Abstract 

The needs to assess robust performances for complex systems and to answer tighter 
regulatory processes (security, safety, environmental control, and health impacts, etc.) 
have led to the emergence of a new industrial simulation challenge: to take uncertainties 
into account when dealing with complex numerical simulation frameworks. Therefore, 
a generic methodology has emerged from the joint effort of several industrial companies 
and academic institutions. EDF R&D, Airbus Group and Phimeca Engineering started 
a collaboration at the beginning of 2005, joined by IMACS in 2014, for the develop¬ 
ment of an Open Source software platform dedicated to uncertainty propagation by 
probabilistic methods, named OpenTURNS for Open source Treatment of Uncertainty, 

Risk ’N Statistics. OpenTURNS addresses the specific industrial challenges attached 
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to uncertainties, which are transparency, genericity, modularity and multi-accessibility. 
This paper focuses on OpenTURNS and presents its main features: OpenTURNS is 
an open source software under the LGPL license, that presents itself as a C++ library 
and a Python TUI, and which works under Linux and Windows environment. All the 
methodological tools are described in the different sections of this paper: uncertainty 
quantification, uncertainty propagation, sensitivity analysis and metamodeling. A sec¬ 
tion also explains the generic wrappers way to link OpenTURNS to any external code. 
The paper illustrates as much as possible the methodological tools on an educational 
example that simulates the height of a river and compares it to the height of a dyke 
that protects industrial facilities. At last, it gives an overview of the main developments 
planned for the next few years. 

Key words: OpenTLTRNS , uncertainty, quantification, propagation, estima¬ 
tion, sensitivity, simulation, probability, statistics, random vectors, multivariate distri¬ 
bution, open source, python module, C++ library, transparency, genericity. 

Introduction 

The needs to assess robust performances for complex systems and to answer tighter reg¬ 
ulatory processes (security, safety, environmental control, and health impacts, etc.) have 
led to the emergence of a new industrial simulation challenge: to take uncertainties into 
account when dealing with complex numerical simulation frameworks. Many attempts 
at treating uncertainty in large industrial applications have involved domain-specific 
approaches or standards: metrology, reliability, differential-based approaches, variance 
decomposition, etc. However, facing the questioning of their certification authorities 
in an increasing number of different domains, these domain-specific approaches are no 
more appropriate. Therefore, a generic methodology has emerged from the joint effort 
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of several industrial companies and academic institutions : [l 2H] reviews these past devel¬ 
opments. The specific industrial challenges attached to the recent uncertainty concerns 
are: 

• transparency: open consensus that can be understood by outside authorities and 
experts, 

• genericity: multi-domain issue that involves various actors along the study, 

• modularity: easy integration of innovations from the open source community, 

• multi-accessibility: different levels of use (simple computation, detailed quan¬ 
titative results, and deep graphical analyses) and different types of end-users 
(graphical interface, Python interpreter, and C++ sources), 

• industrial computing capabilities: to secure the challenging number of simula¬ 
tions required by uncertainty treatment. 

As no software was fully answering the challenges mentioned above, EDF R&D, 
Airbus Group and Phimeca Engineering started a collaboration at the beginning of 
2005, joined by EMACS in 2014, for the development of an open Source software 
platform dedicated to uncertainty propagation by probabilistic methods, named Open- 
TURNS for open source Treatment of Uncertainty, Risk ’N Statistics ram Open- 
TURNS is actively supported by its core team of four industrial partners (IMACS 
joined the consortium in 2014) and its industrial and academic users community that 
meet through the web site www. openturns .org and annually during the OpenTURNS 
Users day. At EDF, OpenTURNS is the repository of all scientific developments on this 
subject, to ensure their dissemination within the several business units of the company. 
The software has also been distributed for several years via the integrating platform 
Salome (27J. 
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Presentation of OpenTURNS 

OpenTURNS is an open source software under the LGPL license, that presents itself 
as a C++ library and a Python TUI, and which works under Linux and Windows 
environment, with the following key features: 

• open source initiative to secure the transparency of the approach, and its open¬ 
ness to ongoing Research and Development (R&D) and expert challenging, 

• generic to the physical or industrial domains for treating of multi-physical prob¬ 
lems; 

• structured in a practitioner-guidance methodological approach, 

• with advanced industrial computing capabilities, enabling the use of massive dis¬ 
tribution and high performance computing, various engineering environments, 
large data models etc., 

• includes the largest variety of qualified algorithms in order to manage uncer¬ 
tainties in several situations, 

• contains complete documentation (Reference Guide, Use Cases Guide, User 
manual, Examples Guide, and Developers’ Guide). 

All the methodological tools are described after this introduction in the different sec¬ 
tions of this paper: uncertainty quantification, uncertainty propagation, sensitivity 
analysis and metamodcling. Before the conclusion, a section also explains the generic 
wrappers way to link OpenTURNS to any external code. 

OpenTURNS can be downloaded from its dedicated website www.openturns.org 
which offers different pre-compilcd packages specific to several Windows and Linux 
environments. It is also possible to download the source hies from the SourceForge 
server ( www.sourceforge.net ) and to compile them within another environment: the 
OpenTURNS Developer’s Guide provides advice to help compiling the source hies. 
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At last, OpenTURNS has been integrated for more than 5 years in the major Linux 

distributions (for example debian, ubuntu, redhat and suze). 

The uncertainty management methodology 

The uncertainty management generic methodology [22] is schematized in Figure [lj It 

consists of the following steps: 

• Step A: specify the random inputs X , the deterministic inputs d, the model 
G (analytical, complex computer code or experimental process), the variable of 
interest (model output) Y and the quantity of interest on the output (central 
dispersion, its distribution, probability to exceed a threshold, ...). The funda¬ 
mental relation writes: 

Y = G(X,d)=G(X), (1) 

with X — (X 1 ,...,X d ). 

• Step B: quantify the sources of uncertainty. This step consists in modeling the 
joint probability density function (pdf) of the random input vector by direct 
methods (e.g. statistical fitting, expert judgment) [ I5j . 

• Step B’: quantify the sources of uncertainty by indirect methods using some real 
observations of the model outputs [32] • The calibration process aims to estimate 
the values or the pdf of the inputs while the validation process aims to model 
the bias between the model and the real system. 

• Step C: propagate uncertainties to estimate the quantity of interest. With re¬ 
spect to this quantity, the computational resources and the CPU time cost of 
a single model run, various methods will be applied: analytical formula, ge¬ 
ometrical approximations, Monte Carlo sampling strategies, metamodcl-based 
techniques, ... EH. EH- 
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• Step C’: analyze the sensitivity of the quantity of interest to the inputs in order 


to rank uncertainty sources ra. Pi- 


For each of these steps, OpenTURNS offers a large number of different methods whose 
applicability depend on the specificity of the problem (dimension of inputs, model 
complexity, CPU time cost for a model run, quantity of interest, etc.). 


Step B : Uncertainty 
Quantification 

Modeling with 
probability distribution : 
direct methods, 
statistics, expertise. 

A_i. 





Step C : Uncertainty Propagation 
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Step B': Quantification 
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Fig. 1 . The uncertainty management methodology. 


Main originality of OpenTURNS 

OpenTURNS is innovative in several aspects. Its input data model is based on the 
multivariate cumulative distribution function (CDF). This enables the usual sampling 
approach, as would be appropriate for statistical manipulation of large data sets, but 
also facilitates analytical approaches. Distributions are classified (continuous, discrete, 
elliptic, etc.) in order to take the best benefit of their properties in algorithms. If 
possible, the exact final cumulative density function is determined (thanks to charac- 
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teristic functions implemented for each distribution, the Poisson summation formula, 
the Cauchy integral formula, etc.). 

OpenTURNS explicitly models the dependence with copulas, using the Sklar 
theorem. Furthermore, different sophisticated analytical treatments may be explored: 
aggregation of copulas, composition of functions from R n into R d , extraction of copula 
and marginals from any distribution. 

OpenTURNS defines a domain specific oriented object language for probability 
modelling and uncertainty management. This way, the objects correspond to mathe¬ 
matical concepts and their inter-relations map the relations between these mathemati¬ 
cal concepts. Each object proposes sophisticated treatments in a very simple interface. 

OpenTURNS implements up-to-date and efficient sampling algorithms (Mersenne- 
Twister algorithm, Ziggurat method, the Sequential Rejection Method, etc.). Exact 
Kolmogorov statistics are evaluated with the Marsaglia Method and the Non Central 
Student and Non Central y 2 distribution with the Benton & Krishnamoorthy method. 

OpenTURNS is the repository of recent results of PhD research carried out 
at EDF R&D: for instance the sparse Polynomial Chaos Expansion method based on 
the LARS method )T], the Adaptive Directional Stratification method |25j which is an 
accelerated Monte Carlo sampling technique, and the maximum entropy order-statistics 
copulas m- 

The flooding model 

Throughout this paper, the discussion is illustrated with a simple application model 
that simulates the height of a river and compares it to the height of a dyke that protects 
industrial facilities as illustrated in Figure [2] When the river height exceeds that of the 
dyke, flooding occurs. This academic model is used as a pedagogical example in H- 
The model is based on a crude simplification of the ID hydro-dynamical equations 


of SaintVenant under the assumptions of uniform and constant flowrate and large 


rectangular sections. It consists of an equation that involves the characteristics of the 
river stretch: 


0.6 



( 2 ) 


where the output variable H is the maximal annual height of the river, B is the river 
width and L is the length of the river stretch. The four random input variables Q, K s , 
Z v and Z m are defined in Table [l] with their probability distribution. The randomness 
of these variables is due to their spatio-temporal variability, our ignorance of their true 
value or some inaccuracies of their estimation. 



Fig. 2. The flood example: simplified model of a river. 


Input Description Unit Probability distribution 


Q Maximal annual flowrate m 3 /s Gumbel <5(1.8e 3 ,1014) 


K. 


Strickler coefficient 


Normal 7V(30, 7.5) 


Z v River downstream level m Triangular T(47.6, 50.5,52.4) 


Z m River upstream level m Triangular T(52.5,54.9,57.7) 


Table 1. Input variables of the flood model and their probability distributions. 
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Uncertainty quantification 

Modelling of a random vector 

OpenTURNS implements more than 40 parametric distributions which are continuous 
(more than 30 families) and discrete (more than 10 families), with several sets of pa¬ 
rameters for each one. Some are multivariate, such as the Student distribution or the 
Normal one. 

Moreover, OpenTURNS enables the building of a wide variety of multivariate distribu¬ 
tions thanks to the combination of the univariate margins and a dependence structure, 
the copula, according to the Sklar theorem: F(xi ,..., Xj) = C (Fi(xi ),..., Fd(xd)) 
where F t is the CDF of the margin X t and C : [0, l] d —> [0,1] the copula. 

OpenTURNS proposes more than 10 parametric families of copula: Clayton, Frank, 
Gumbel, Farlie-Morgenstein, etc. These copula can be aggregated to build the copula 
of a random vector whose components are dependent by blocks. Using the inverse re¬ 
lation of the Sklar theorem, OpenTURNS can extract the copula of any multivariate 
distribution, whatever the way it has been set up: for example, from a multivariate 
distribution estimated from a sample with the kernel smoothing technique. 

All the distributions can be truncated in their lower and/or upper area. In addition to 
these models, OpenTURNS proposes other specific constructions. Among them, note 
the random vector which writes as a linear combination of a finite set of independent 
variables: X = ao + d\X\ + .. .a^X n thanks to the python command, written for 
N — 2 with explicit notations: 

»>myX= RandomMixture ([ distXl , distX2] , [al , a2] , aO) 

In that case, the distribution of X is exactly determined, using the characteristic 
functions of the X, distributions and the Poisson summation formula. 

OpenTURNS also easily models the random vector whose probability density func¬ 
tion (pdf) is a linear combination of a finite set of independent pdf: fj£ = a ifx i + 
• • • a N fx N thanks to the python command, with the same notations as previously (the 
weights are automatically normalized): 


»>mypdfX= Mixture ([ distXl , distX2 ] , [al , a2 ]) 
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Moreover, OpenTURNS implements a random vector that writes as the random 
sum of univariate independent and identically distributed variables, this randomness 
being distributed according to a Poisson distribution: X = X^Xj, N r s./ P(A), 
thanks to the python command: 

»>d= CompoundDistribution (lambda, distX) 

where all the variables X, are identically distributed according to distX. In that 
case, the distribution of X is exactly determined, using the characteristic functions of 
the Xi distributions and the Poisson summation formula. 

In the univariate case, OpenTURNS exactly determines the pushforward distribution 
T> of any distribution T >o through the function / : M. —y M, thanks to the python 
command (with straight notations): 

»>d= CompositeDistribution(f , dO) 

Finally, OpenTURNS enables the modeling of a random vector (Xi,...,X d ) 
which almost surely verifies the constraint X = X\ < ■ ■ ■ < X d , proposing a copula 
adapted to the ordering constraint [S|. OpenTURNS verifies the compatibility of the 
margins F t with respect to the ordering constraint and then builds the associated 
distribution, thanks to the python command, written in dimension 2: 

»>d=MaximumEntropyOrderStatisticsDistribution ([ distXl , distX2 ]) 

Figure [3] illustrates the copula of such a distribution, built as the ordinal sum 

of some maximum entropy order statistics copulae. 

The OpenTURNS python script to model the input random vector of the tutorial 
presented previously is as follows: 

#Margin distributions : 

»>dist_Q = Gumbel (1.8 e —3, 1014) 

»>dist_Q = TruncatedDistribution (dist_Q , 0.0 , TruncatedDistribution .LOWER) 
»>dist_K = Normal (30.0, 7.5) 

»>dist_K = TruncatedDistribution (dist_K , 0 . , TruncatedDistribution .LOWER) 
»>dist_Zv = Triangular (47.6 ,50.5 ,52.4) 

»>dist__Zm = Triangular ( 5 2.5 ,5 4.9 ,5 7.7) 
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Fig. 3. An example of maximum entropy copula which almost surely satisfies the ordering constraint: 


Xi < X 2 . 


# Copula in dimension 4 f or (Q,K,Zv,Zm) 

»>R=C orrelationMatrix(2) 

»>R[0,1] =0.7 

»>copula = ComposedCopula ([ IndependentCopula (2) , NormalCopula (R) ]) 

# Final distribution for (Q,K, Zv ,Zm) 

»>distInput=ComposedDistribution ([loi_Q , loi_K , loi_Zv , loi_Zm], copula) 

# Final random vector (Q,K, Zv ,Zm) 

»>inputVector=RandomVector (clistlnput) 

Note that OpenTURNS can truncate any distribution to a lower, an upper 
bound or a given interval. Furthermore, a normal copula models the dependence be¬ 
tween the variables Z v and Z m , with a correlation of 0.7. The variables ( Q,K ) are 
independent. Both blocks ( Q,K) and (Z v ,Z m ) are independent. 

Stochastic processes 

OpenTURNS implements some multivariate random fields X : i? x T> —> where 

D G R s is discretized on a mesh. The User can easily build and simulate a random 
walk, a white noise as illustrated in Figures [4] and [5} The python commands write: 

»XnyWN = WhiteNoise (myDist , myMesh) 

»XnyRW = RandomWalkf my Origin , myDist, myTimeGrid) 









12 


2D Random Walk with U({-1,1}) steps 


o.o 



0.5 


Field with linear interpolation 


1.0 


1.5 


2.0 


X 



-40 


-20 


20 


40 


Fig. 4. A Normal bivariate white noise. 


XI 


Fig. 5. A Normal bivariate random walk. 


Any field can be exported into the VTK format which allows it to be visualized 
using e.g.ParaView (www.paraview.org). 

Multivariate ARMA stochastic processes X : Q x [0, T] —>■ are implemented in 

OpenTURNS which enables some manipulations on times series such as the Box Cox 
transformation or the addition / removal of a trend. Note that the parameters of the 
Box Cox transformation can be estimated from given fields of the process. 
OpenTURNS models normal processes, whose covariance function is a parametric 
model (e.g. the multivariate Exponential model) as well as defined by the User as 
illustrated in Figure [6j Stationary processes can be defined by its spectral density 
function (e.g. the Cauchy model). 


C as a function of (s,t) 



-4 


-2 


Fig. 6. A User Defined non stationary covariance function and its estimation from several given fields. 
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With explicit notations, the following python commands create a stationary 
Normal process defined by its covariance function, discretized on a mesh, with an 
additional trend: 

»>myNormalProcess=TemporalNormalProcess (myTrend , myCovarianceModel , myMesh) 

Note that OpenTURNS enables the mapping of any stochastic processes X into 
a process Y through a function f: Y — f(X) where the function / can consist, for 
example, of adding or removing a trend, applying a Box Cox transformation in order 
to stabilize the variance of X. The python command is, with explicit notations: 

»>myYprocess=CompositeProcess ( f , myXprocess) 

Finally, OpenTURNS implements multivariate processes defined as a linear com¬ 
bination of K deterministic functions '■ IR dl lU 2 : 

I< 

X(u,x) = 

i =1 

where {A \,..., A K ) is a random vector of dimension K. The python command writes: 
»>myX =FunctionalBasisProcess (myRandomCoeff, myBasis , myMesh) 


Statistics estimation 

OpenTURNS enables the User to estimate a model from data, in the univariate as 
well as in the multivariate framework, using the maximum likelihood principle or the 
moments based estimation. 

Some tests, such as the Kolmogorov-Smirnov test, the Chi Square test and the Anderson 
Darling test (for normal distributions), are implemented and can help to select a model 
amongst others, from a sample of data. The python command to build a model and 
test it, writes: 

»>estimatedBeta = BetaFactory (sample) 

»>testResult = FittingTest . Kolmogorov (sample , estimatedBeta) 

OpenTURNS also implements the kernel smoothing technique which is a non- 
parametric technique to fit a model to data: any distribution can be used as kernel. 
In the multivariate case, OpenTURNS uses the product kernel. It also implements 
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an optimized strategy to select the bandwidth, depending on the number of data in 
the sample, which is a mix between the Silverman rule and the plugin one. Note that 
OpenTURNS proposes a special treatment when the data are bounded, thanks to the 
mirroring technique. The python command to build the non-parametric model and to 
draw its pdf, is as simple as the following one: 

»>estimatedDist = KernelSmoothing (). build (sample ) 

»>pdfGraph = estimatedDist .drawPDF() 

Figure [7] illustrates the resulting estimated distributions from a sample of size 
500 issued from a Beta distribution: the Kernel Smoothing method takes into account 
the fact that data are bounded by 0 and 1. The histogram of the data is drawn to 
enable comparison. 

Several visual tests are also implemented to help to select models: among them, 
the QQ-plot test and the Henry line test which writes (in the case of a Beta distributino 
for the QQ-plot test): 

»>graphQQplot = VisualTest .DrawQQplot(sample , Beta()) 

»>graphHenryLine = VisualTest . DrawHenryLine (sample ) 

Figure [8] illustrates the QQ-plot test on a sample of size 500 issued from a Beta 

distribution: the adequation seems satisfying! 

Stochastic processes also have estimation procedures from sample of fields or, 
if the ergodic hypothesis is verified, from just one held. Multivariate ARMA processes 
are estimated according to the BIC and AIC criteria and the Whittle estimator, which 
is based on the maximization of the likelihood function in the frequency domain. The 
python command to estimate an ARMA(p, q) process of dimension d, based on a 
sample of time series, writes: 

»>estimatedARMA = ARMALikelihood (p , q , d ). build ( sampleTimeSeries ) 

Moreover, OpenTURNS can estimate the covariance function and the spectral 
density function of normal processes from given fields. For example, the python com- 


15 


Sample fitting 



Sample versus Beta: QQPIot test 



Fig. 7. Beta distribution estimation from a sam¬ 
ple of size 500: parametric estimation versus ker¬ 
nel smoothing technique. 


Fig. 8. QQ-plot test: theoretical model Beta ver¬ 
sus the sample of size 500. 


mand to estimate a stationary covariance model from a sample of realizations of the 
process: 

»>myCovFunc = S tat ionary C ovar ianceModelFac t ory (). build (sampleProcess) 


This estimation is illustrated in Figure [6} 


Conditioned distributions 

OpenTURNS enables the modeling of multivariate distributions by conditioning. Sev¬ 
eral types of conditioning are implemented. 

At first, OpenTURNS enables the creation of a random vector X whose distribution 
T^X\0 w h° se parameters 0 form a random vector distributed according to the distri¬ 
bution Dq . The python command writes: 

»>myXrandVect = ConditionalRandomVector ( distXgivenTheta , distTheta) 

Figure [9] illustrates a random variable X distributed according to a Normal 
distribution: T^X\0=(m e) = A Iormal(M, E), which parameters are defined by M ~ 
Uniform([ 0,1]) and X ~ Exponential (A = 4). The probability density function of A" 
has been fitted with the kernel smoothing technique from n = 10 6 realizations of X 
with the normal kernel. It also draws, for comparison needs, the probability density 
function of X in the case where the parameters are fixed to their mean value. 
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Furthermore, when the random vector 0 is defined as © = g(Y) where the random 
vector follows a known distribution T>y and g is a given function, OpenTURNS creates 
the distribution of X with the python command: 

»>finalDist = ConditionalDistribution (distXGgivenTheta , distY , g) 

Figure [To| illustrates the distribution of A" that follows a Uniform(A, B) dis¬ 
tribution, with (A, B) = g(Y ), g : R —» R 2 , g(Y) = (Y, 1 + Y 2 ) and Y follows a 
Uniform^— 1,1) distribution. 


Distribution Y 



Fig. 9. Normal distribution with random or fixed 
parameters. 


Conditional distribution 



x 


Fig. 10. UniformiY, 1 + Y 2 ), with Y ~ 
Uniform (—1,1). 


Bayesian Calibration 

Finally, OpenTURNS enables the calibration of a model (which can be a computer 
code) thanks to the Bayesian estimation, which is the evaluation of the model’s pa¬ 
rameters. More formally, let’s consider a model G that writes: y = G(x,0) where 
x £ R dl , y £ R ds and 6 £ R d2 is the vector of unknown parameters to calibrate. 
The Bayesian calibration consists in estimating 0, based on a certain set of n inputs 
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(a? 1 , ..., x n ) (an experimental design) and some associated observations (z 1 , ..., z n ) 
which are regarded as the realizations of some random vectors (Z 1 ,..., Z n ), such that, 
for all i, the distribution of Z l depends on y' = g(x l ,6). Typically, Z l = Y l + e l 
where e % is a random measurement error. Once the User has defined the prior distribu¬ 
tion of 6, OpenTURNS maximizes the likelihood of the observations and determines 
the posterior distribution of 6 , given the observations, using the Metropolis-Hastings 
algorithm jHJESj- 


Uncertainty propagation 

Once the input multivariate distribution has been satisfactorily chosen, these uncer¬ 
tainties can be propagated through the G model to the output vector Y. Depending on 
the final goal of the study (min-max approach, central tendency, or reliability), several 
methods can be used to estimate the corresponding quantity of interest, tending to 
respect the best compromise between the accuracy of the estimator and the number of 
calls to the numerical, and potentially costly, model. 

Min-Max approach 

The aim here is to determine the extreme (minimum and maximum) values of the 
components of Y for the set of all possible values of X. Several techniques enable it to 
be done : 

• techniques based on design of experiments : the extreme values of Y are sought 
for only a finite set of combinations (x 1 ,..., x”), 

• techniques using optimization algorithms. 

Techniques based on design of experiments 

In this case, the min-max approach consists of three steps: 
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• choice of experiment design used to determine the combinations (x 1 ,..., x") of 
the input random variables, 

• evaluation of y* = G(x l ) for i = 1,..., N, 

• evaluation of min i<i<NVi and of maxjxj <NVi, together with the combina¬ 
tions related to these extreme values: Xfc imin = argmin 1<i<jV y^ and Xfc ]max = 
argmax^^yf. 

The type of design of experiments impacts the quality of the rneta model and 
then on the evaluation of its extreme values. OpenTURNS gives access to two usual 
family of design of experiments for a min-max study : 

• some stratified patterns (axial, composite, factorial or box patterns) Here are 
the two command lines that generate a sample from a 2-level factorial pattern. 

»>myCenteredReductedGrid = Factorial (2, levels) 

»>mySample = myCenteredReducedGrid . generate () 

• some weighted patterns that include on the one hand, random patterns (Monte 
Carlo, LHS), and on the other hand, low discrepancy sequences (Sobol, Faure, 
Halton, Reverse Halton and Haselgrove, in dimension n > 1). The following 
lines illustrates the creation of a Faure sequence in diension 2 or a Monte Carlo 
design experiments form a bi-dimensional Normal(0,l) distribution. 

# Sobol Sequence Sampling 

»>mySobolSample = FaureSequence ( 2 ). generate (1000) 

# Monte Carlo Sampling 

»>myMCSample = MonteCarloExperiment (Normal (2) , 100) 

Figures [TT] and [12] respectively illustrate a design of experiments issued from a 
Faure sequence or a Normal distribution in dimension 2. 

Techniques based on optimization algorithm 
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Sequence of 1000 points 




Fig. 11. The first 1000 points according of a Fig. 12. Sample of 1000 points according of a 
Faure sequence of dimension 2. Normal(0,l) distribution of dimension 2. 

In this kind of approach, the min or max value of the output variable is sought 
thanks to an optimization algorithm. OpenTURNS offers several optimization algo¬ 
rithms for the several steps of the global methodology. Here the TNC (Truncated 
Newton Constrainted) is often used, which minimizes a function with variables subject 
to bounds, using gradient information. More details may be found in j2SJ- 

# For the research of the min value 

»>myAlgoTNC = TNC( TNCSpecificParameters () , limitStateFunction , 

intervalOptim , startingPoint , TNC.MINIMIZATION) 

# For the research of the max value 

»>myAlgoTNC = TNC( TNCSpecificParameters () , limitStateFunction , 

intervalOptim , startingPoint , TNC.MAXIMIZATION) 

# Run the research and extract the results 
»>myAlgoTNC. run () 

»>myAlgoTNCResult = BoundConstrainedAlgorithm (myAlgoTNC). getResult () 
»>optimalValue = myAlgoTNCResult. getOptimalValue () 


Central tendency 

A central tendency evaluation aims at evaluating a reference value for the variable 
of interest, here the water level H , and an indicator of the dispersion of the variable 
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around the reference. To address this problem, mean /iy = e(Y'), and the standard 
deviation cry = \Jy(Y) of Y are here evaluated using two different methods. 

First, following the usual method within the Measurement Science community 
ra, /iy and cry have been computed under a Taylor first order approximation of the 
function Y = G(X) (notice that the explicit dependence on the deterministic variable 
d is here omitted for simplifying notations): 


//y~G(E(X)) 
A A dG 

1=1 J = 1 L 


dG 


<X)dXj 


*X) 


Pij&i&j i 


(3) 

(4) 


cjj and crj being the standard deviation of the ith and jth component X t and Xj 
of the vector X and p^ their correlation coefficient. Thanks to the formulas above, the 
mean and the standard deviation of H are evaluated as 52.75m and 1.15m respectively. 


»>myQuadCum = QuadraticCumul ( output Variable ) 

# First order Mean 

»>meanFirstOrder = myQuadCum. getMeanFirstOrder () [0 ] 

# Second order Mean 

»>meanSecondOrder = myQuadCum. getMeanSecondOrder ()[0] 

# First order Variance 

»>varFirst Order = myQuadCum. getCovariance ()[ 0 , 0 ] 


Then, the same quantities have been evaluated by a Monte Carlo evaluation : a 
set of 10000 samples of the vector X is generated and the function G(X) is evaluated, 
thus giving a sample of H. The empirical mean and standard deviation of this sample 


are 52.75m and 1.42m respectively. Figure 13 shows the empirical histogram of the 
generated sample of H. 


# Create a random sample of the output variable of interest of size 10000 
»>outputSample = output Variable . getNumericalSample (10000) 

# Get the empirical mean 

»>empiricalMean = outputSample . computeMean () 
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sample histogram 



Fig. 13. Empirical histogram of 10000 samples of H. 

# Get the empirical covariance matrix 

»>empiricalCovarianceMatrix = outputSample . computeCovariance () 


Failure probability estimation 

This section focuses on the estimation of the probability for the output Y to exceed a 
certain threshold s, noted Pf in the following. If s is the altitude of a flood protection 
dyke, then the above excess probability, Pf can be interpreted as the probability of an 
overflow of the dyke, i.e. a failure probability. 

Note that an equivalent way of formulating this reliability problem would be 
to estimate the (1 — p)-th quantile of the output’s distribution. This quantile can 
be interpreted as the flood height q p which is attained with probability p each year. 
T = 1/p is then seen to be a return period, i.e. a flood as high than qi/p occurs on 
average every T years. 

Hence, the probability of overflowing a dyke with height s is less than p (where 
p, for instance, could be set according to safety regulations) if and only if s > q p , i.e. 
if the dyke’s altitude is higher than the flood with return period equal to T = 1/p. 
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FORM 

A way to evaluate such failure probabilities is through the so-called First Order Relia¬ 
bility Method (FORM) [9]. This approach allows, by using an equiprobabilistic trans¬ 
formation and an approximation of the limit-state function, the evaluation with a much 
reduced number of model evaluations, of some low probability as required in the reliabil¬ 
ity held. Note that OpenTURNS implements the Nataf transformation where the input 
vector X has a normal copula, the generalized Nataf transformation when X has an 
elliptical copula, and the Rosenblatt transformation for any other cases dEi nsnaim. 

The probability that the yearly maximal water height H exceeds s=58m is 
evaluated using FORM. The Hasofer-Lind Reliability index was found to be equal to: 
(3hl = 3.04, yielding a final estimate of: 


P 


f,FORM 


= 1.19 x 10" 3 . 


The method gives also some importance factors that measure the weight of each 


input variable in the probability of exceedance, as shown on Figure 14 


FORM Importance Factors - Event Zc > 58.0 



Fig. 14. FORM Importance Factors. 


»>myFORM = PORM( Cobyla () , myEvent, meanlnputVector) 
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»>myFORM. run () 

»>FormResult = myFORM. getResult () 

»>pFORM = FormResult. getEventProbability () 

»>HasoferIndex = FormResult . getHasoferReliabilitylndex () 

# Importance factors 

»>importanceFactorsGraph = FormResult. drawImportanceFactors () 


Monte Carlo 


Whereas the FORM approximation relies on strong assumptions, the Monte Carlo 
method is always valid, independently from the regularity of the model. It is neverthe¬ 
less much more computationally intensive, covering all the input domain to evaluate 
the probability of exceeding a threshold. It consists in sampling many input values 
from the input vector joint distribution, then computing the correspond¬ 
ing output values = g(X^). The excess probability Pf is then estimated by the 
proportion of sampled values that exceed t : 

1 N 

p f,MC = jr 1 {y( i )> s }- ( 5 ) 

i= 1 

The sample average of the estimation error Pf.MC — Pf decreases as 1 / x/iV, and can 
be precisely quantified by a confidence interval derived from the central limit theorem. 
In the present case the results are: 


P fM c = 1-50 x 1(T 3 , 


with the following 95% confidence interval: 


l P f ,MC 


1.20 x 10 -3 ,1.79 x 10 


-3 


These results are coherent with those of the FORM approximation, confirming that 
the assumptions underlying the latter are correct. Figure [15] shows the convergence of 
the estimate depending on the size of the sample, obtained with OpenTURNS . 
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MonteCarlo convergence graph at level 0.95 



Fig. 15. Monte Carlo Convergence graph. 

»>myEvent = Event ( output Variable , GreaterQ, threshold) 

»>myMonteCarlo = MonteCarlo (myEvent) 

# Specify the maximum number of simulations 
>»numberMaxSimulation = 100000 

»>myMonteCarlo . setMaximumOuterSampling (numberMaxSimulation ) 

# Perform the algorithm 
»>myMonteCarlo. run () 

# Get the convergence graph 

»>convergenceGraph = myMonteCarlo . drawProbabilityConvergence () 
»>convergenceGraph . draw (" MonteCarloCovergenceGraph ") 

Importance Sampling 

An even more precise estimate can be obtained through importance sampling [51] . 
using the Gaussian distribution with identity covariance matrix and mean equal to the 
design point u* as the proposal distribution. Many values (C^)i<i<jv are sampled from 
this proposal. Because <f n {u — u*) is the proposal density from which the have been 
sampled, the failure probability can be estimated without bias by: 
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p _ i v, *.( C/W ) 

iv i= 1 


( 6 ) 


'<A„(t/w-«•) 

The rationale of this approach is that by sampling in the vicinity of the failure domain 
boundary, a larger proportion of values fall within the failure domain than by sam¬ 


pling around the origin, leading to a better evaluation of the failure probability, and a 
reduction in the estimation variance. Using this approach,the results are : 


P fJ s = 1.40 x 10~ 3 


As in the simple Monte-Carlo approach, a 95%-level confidence interval can be derived 
from the output of the Importance Sampling algorithm. In the present case, this is 
equal to: 


J PfdS 


1.26 x 1CT 3 ,1.53 x 10“ 3 


and indeed provides tighter confidence bounds for Pf. 

# Specify the starting point from FORM algorithm 
»>standardPoint = FormResult . getStandardSpaceDesignPoint () 

# Define the importance distribution 
»>sigma = [1.0, 1.0, 1.0, 1.0] 

»>importanceDistrib = Normal(standardPoint , sigma , CorrelationMatrix (4)) 

# Define the IS algorithm : event, distribution , criteria of convergence 
»>myAlgoImportanceSampling = ImportanceSampling (myStandardEvent , importanceDistrib ) 

»>myAlgoImportanceSampling . setMaximumOuterSampling (maximumOuterSampling_IS) 

»>myAlgoImportanceSampling . setMaximumCoefficientOfVariation (0.05) 


Directional Sampling 

The directional simulation method is an accelerated sampling method, that involves 
as a first step a preliminary iso-probabilistic transformation as in FORM method. The 
basic idea is to explore the space by sampling in several directions in the standard 
space. The final estimate of the probability Pf after N simulations is the following: 
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1 N 

Pf,DS = T7 H 9i 

JV i =1 

where g* is the probability obtained in each explored direction. A central limit theorem 
allows to access to some confidence interval on this estimate. More details on this 
specific method can be found in [32]. 

In practice in OpenTURNS , the Directional Sampling simulation requires the 
choice of several parameters in the methodology : a sampling strategy to choose the 
explored directions, a "root strategy" corresponding to the way to seek the limit state 
function (i.e. a sign change) along the explored direction and a non-linear solver to 
estimate the root. A default setting of these parameters allows the user to test the 
method in one command line : 

»>myAlgo = DirectionalSampling (myEvent) 


Subset Sampling 

The subset sampling is a method for estimating rare event probability, based on the 
idea of replacing rare failure event by a sequence of more frequent events F t . 

Fi D F 2 D ■ ■ • D F m = F 

The original probability is obtained conditionaly to the more frequent events : 

m m 

Pf = p(F m ) = p(n *;) = c(n) n miu-i) 

i =1 i=2 

In practice, the subset simulation shows a substantial improvement (N? rv./ log Pf) 
compared to crude Monte Carlo (Nt ~ p f ) sampling when estimating rare events. 
More details on this specific method can be found in |2J. 

OpenTURNS provides this method through a dedicated module. Here also, some 
parameters of the methods have to be chosen by the user : a few command lines allows 
the algorithm to be set up before its launch. 


>»mySSAlgo=SubsetSampling (myEvent) 
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# Change the target conditional probability of each subset domain 
»>mySSAlgo. setTargetProbability (0.9) 

# Set the width of the MCMC random walk uniform distribution 
»>mySSAlgo. setProposalRange (1-8) 

# This allows to control the number of samples per step 
»>mySSAlgo. setMaximumOuterSampling (10000) 

# Run the algorithm 
»>mySSAlgo . run () 


Sensitivity analysis 

The sensitivity analysis aims to investigate how a given computational model answers 
to variations in its inputs. Such knowledge is useful for determining the degree of resem¬ 
blance of a model and a real system, distinguishing the factors that mostly influence 
the output variability and those that are insignificant, revealing interactions among in¬ 
put parameters and correlations among output variables, etc. A detailed description of 
sensitivity analysis methods can be found in [36] and in the Sensitivity analysis chapter 
of the Springer Handbook. In the global sensitivity analysis strategy, the emphasis is 
put on apportioning the output uncertainty to the uncertainty in the input factors, 
given by their uncertainty ranges and probability distributions. 


Graphical tools 


In sensitivity analysis, graphical techniques have to be used first. With all the scatter- 
plots between each input variable and the model output, one can immediately detect 
some trends in their functional relation. The following instructions allow scatterplots of 


Figure 16 to be obtained from a Monte Carlo sample of size N = 1000 of the flooding 
model. 


»> inputSample = inputRandomVector . getNumericalSample (1000) 
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»> inputSample . setDescription ([ ’Q’ , ’K’, ’Zv’, ’Zm’]) 

»> outputSample = finalModelCrue (inputSample ) 

»> outputSample . setDescription ( [ ’H’ ]) 

# Here, stack both samples in one 
»> inputSample . stack (outputSample ) 

»>myPairs = Pairs (inputSample ) 

»>myGraph = Graph () 

»>myGraph . add (myPairs ) 


In the right column of Figure 16 it is clear that the strong and rather linear 
effects of Q and Z v on the output variable H. In the plot of third line and fourth 
column, it is also clear that the dependence between Z v and Z m , which comes from the 
large correlation coefficient introduced in the probabilistic model. 
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Fig. 16. Scatterplots between the inputs and the output of the flooding model: each combination 
(input i, input j) and (input i, output) is drawn, which enables to exhibit some correlation patterns. 


However scatterplots do not capture some interaction effects between the inputs. 
Cobweb plots are then used to visualize the simulations as a set of trajectories. The 
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following instructions allow the cobweb plots of Figure 17 to be obtained where the 
simulations leading to the largest values of the model output H have been colored in 
red. 


»>inputSample = inputRandomVector . getNumericalSample (1000) 
>»outputSample = finalModelCrue (inputSample) 

# Graph 1 : value based scale to describe the Y range 
»>minValue = outputSample . computeQuantilePerComponent (0.05) [0] 
»>maxValue = outputSample . computeQuantilePerComponent (0.95) [0] 
»>myCobweb = VisualTest . DrawCobWeb( inputSample , outputSample, 

minValue , max Value , ’red’, False) 


The cobweb plot allows us to immediately understand that these simulations 
correspond to large values of the flowrate Q and small values of the Strickler coefficient 
K s . 



Fig. 17. Cobweb plot for the flooding model: each simulation is drawn. The input marginal values 
are linked to the output value (last column). All the simulations that led to a high quantile of the 
output are drawn in red: the cobweb plot enables to detect all the combinations of the inputs they 


come from. 
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Sampling-based methods 

In order to obtain quantitative sensitivity indices rather than qualitative information, 
one may use some sampling-based methods which often suppose that the input variables 
are independent. The section illustrates some of these methods on the flooding model 
with independence between its input variables. 

If the behavior of the output Y compared to the input vector X is overall linear, 
it is possible to obtain quantitative measurements of the inputs influences from the 
regression coefficients cq of the linear regression connecting Y to the X = (Ad, ..., X p ). 
The Standard Regression Coefficient (SRC), defined by 

SRC, = cq— (for i = 1... p), (7) 

cry 

with Oi (resp. ay) the standard deviation of A \ (resp. Y), measures the variation of 
the response for a given variation of the parameter Ad. In practice, the coefficient R 2 
(the variance percentage of the output variable Y explained by the regression model) 
also helps to check the linearity: if R 2 is close to one, the relation connecting Y to all 

the parameters Ad is almost linear and the SRC sensitivity indices make sense. 

The following instructions allow the results of Table [2] to be obtained from a 
Monte Carlo sample of size N = 1000. 

»>inputSample = inputRandomVector . getNumericalSample (1000) 

>»outputSample = finalModelCrue (inputSample ) 

»>SRCCoefficient = CorrelationAnalysis . SRC( inputSample , outputSample ) 
»>linRegModel=LinearModelFactory () . build (inputSample , outputSample ,0.90) 
»>Rsquared = LinearModelTest . LinearModelRSquared (inputSample , 

outputSample , linRegModel ,0.90) 

The SRC values confirm our first conclusions drawn from the scatterplots vi¬ 
sual analysis. As R 2 = 0.97 is very close to one, the model is quasi-linear. The SRC 


coefficients are sufficient to perform a global sensitivity analysis. 
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Table 2. Regression coefficients and SRC of the flood model inputs (ao = —0.1675 and Ft? = 0.97). 



Q 

K s 

Z v 

Z m 

Oii 

3.2640 

0.0012 

-0.0556 

1.1720 

SRC; 

0.3462 

0.0851 

0.6814 

0.0149 


Several other estimation methods are available in OpenTURNS for a sensitivity 
analysis purpose: 

• derivatives and Pearson correlation coefficients (linearity hypothesis between 
output and inputs), 

• Spearman correlation coefficients and Standard Rank Regression Coefficients 
(monotonicity hypothesis between output and inputs), 

• reliability importance factors with the FORM/SORM importance measures pre¬ 
sented previously (section ), 

• variance-based sensitivity indices (no hypothesis on the model). These last in¬ 
dices, often known as Sobol’ indices and defined by 

Si = — -y -—j- — (first order index) and Sr, = ^ Sj+ y~] S tJ +... (total index), 
V ar O ) i =1 i<j 

( 8 ) 

are estimated in OpenTURNS with the classic pick-freeze method based on two 
independent Monte Carlo samples pH]- In OpenTURNS, other ways to compute 
the Sobol’ indices are the Extended FAST method [31] and the coefficients of 
the polynomial chaos expansion [38]. 

Metamodels 

When each model evaluation is time consuming, it is usual to build a surrogate model 
which is a good approximation of the initial model and which can be evaluated at 


negligible cost. OpenTURNS proposes some usual polynomial approximations: the lin- 
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ear or quadratic Taylor approximations of the model at a particular point, or a linear 
approximation based on the least squares method. Two recent techniques implemented 
in OpenTURNS are detailed here: the polynomial chaos expansion and the kriging 
approximation. 

Polynomial chaos expansion 

The polynomial chaos expansion enables the approximation of the output random 
variable of interest Y = G(X ) with g : by the surrogate model: 

Y = E » k*k ° T(X) 

k&K 

where otk G M p , T is an isoprobabilistic transformation (e.g. the Rosenblatt transforma¬ 
tion) which maps the multivariate distribution of X into the multivariate distribution 
H = nti/h, an d (^fc)fceN a multivariate polynomial basis of Cj l (R d ,W > ) which is or¬ 
thonormal according to the distribution ft. K is a finite subset of N. Y is supposed to 
be of finite second moment. 

OpenTURNS proposes the building of the multivariate orthonornal basis (^(cc))*,^ 
as the cartesian product of orthonormal univariate polynomial family : 

M z ) = Vfriz 1 ) * 2 ) * ''' * &k d (zd) 

The possible univariate polynomial families associated to continuous measures 

are : 

• Hcrmite, which is orthonormal with respect to the Norma,l( 0,1) distribution, 

• Jacobi(a, /3, pararn ), which is orthonormal with respect to the Beta{(3 + 1, a + 
P + 2, —1,1) distribution if pararn = 0 (default value) or to the Beta(a, P, —1,1) 
distribution if pararn = 1, 

• Laguerre(fc), which is orthonormal with respect to the Gamma(k + 1,1, 0) dis¬ 


tribution, 
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• Legendre, which is orthonormal with respect to the Uniform(—l, 1) distribu¬ 
tion. 


OpenTURNS proposes three strategies to truncate the multivariate orthonormal 
basis to the finite set K : these strategies select different terms from the multivariate 
basis, based on a convergence criterion of the surrogate model and the cleaning of the 
less significant coefficients. 

The coefficients of the polynomial decomposition writes: 


ol = argmin aeR KE M 


goT-\Z)- Y,<x k MZ) 

k&K 


( 9 ) 


as well as: 


cx=[E,\goT-\Z)E k {Z)]) k (10) 

where Z = T(X ) is distributed according to /i. 

It corresponds to two points of view implemented by OpenTURNS : the relation (|9]) 
means that the coefficients ( ak)keK minimize the mean quadratic error between the 


model and the polynomial approximation; the relation (10) means that a is the scalar 
product of the model with the k — th element of the orthonormal basis (\Pk)keK- In 
both cases, the expectation is approximated by a linear quadrature formula that 
writes, in the general case: 


EAS(z)}~Y.^S(Si) (ii) 

iei 

where / is a function Li(/i). The set /, the points (Sj) ie j and the weights (u l ) ie [ are 
evaluated from weighted designs of experiments which can be random (Monte Carlo 
experiments, and Importance sampling experiments) or deterministic (low discrepancy 
experiments, User given experiments, and Gauss product experiments). 

At last, OpenTURNS gives access to: 


• the composed model h : Z i—>■ Y = G o T 1 (Z), which is the model of the 

reduced variables Z. Then h = y, oLk&k, 

km 
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# the coefficients of the polynomial approximation : (ctk)k&K, 

# the composed meta model: h, which is the model of the reduced variables reduced 
to the truncated multivariate basis (&k)k£K- Then h = ^ cxk&k, 

k&K 

# the meta model: g : X i —y Y — h o T(X) which is the polynomial chaos approx¬ 
imation as a NumericalMathFunction. Then g = ° T - 

k&K 

When the model is very expensive to evaluate, it is necessary to optimize the 
number of coefficients of the polynomials chaos expansion to be calculated. Some spe¬ 
cific strategies have been proposed by [6] for enumerating the infinite polynomial chaos 
series: OpenTURNS implements the hyperbolic enumeration strategy which is inspired 
by the so-called sparsity-of-effects principle. This strategy states that most models 
are principally governed by main effects and low-order interactions. This enumeration 

strategy selects the multi-indices related to main effects. 

The following lines illustrates the case where the model Q : x H» xsinrc and 
where the input random vector follows a Uniform distribution on [—1,1] 

# Define the model 

»> model = NumericalMathFunction ([ ’x ’] , [ ’x* sin (y) ’ ]) 

# Create the input distribution 
»> distribution = Uniform () 

# Construction of the orthonormal basis 
»> polyColl = [0.] 

»> polyColl [0] = StandardDistributionPolynomialFactory ( distribution ) 

»> enumerateFunction = LinearEnumerateFunction (1) 

»> productBasis = OrthogonalProductPolynomialFactory ( poly Coll , enumerateFunction) 

# Truncature strategy of the multivariate orthonormal basis 

# Choose all the polynomials of degree <= 4 
»> degree = 4 

»> indexMax = enumerateFunction . getStrataCumulatedCardinal ( degree ) 

# Keep all the p olynomials of degree <= 4 
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# which corresponds to the 5 first ones 

»> adaptiveStrategy = FixedStrategy (productBasis , indexMax) 

# Evaluation strategy of the approximation c o efficients 
»> samplingSize = 50 

»> experiment = MonteCarloExperiment ( samplingSize ) 

»> projectionStrategy = LeastSquaresStrategy ( experiment) 

# Creation of the Functional Chaos Algorithm 

»> algo = FunctionalChaos Algorithm (model, distribution , adaptiveStrategy, 
... projectionStrategy) 

»> algo, run () 

# Get the result 

»> functionalChaosResult = algo . getResult () 

»> metamodel = functionalChaosResult . getMetaModel () 


Figure 18 illustrates the result. 


Polynomial Chaos Expansion 



X 


Fig. 18. An example ofa polynomial chaos expansion: the blue line is the reference function Q : x <—> 
a; sin a; and the red one its approximation only valid on [—1,1]. 


The Kriging approximation 


Kriging (also known as Gaussian process regression) [33, [371 3D, 23] is a Bayesian 
technique that aims at approximating functions (most often in order to surrogate them 
because they are expensive to evaluate). In the following it is assumed the aim is to 
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surrogate a scalar-valued model G : x i-» y. Note the OpenTURNS implementation of 
Kriging can deal with vector-valued functions (G : x i —> y), with simple loops over each 
output. It is also assumed the model was run over a design of experiments in order to 
produce a set of observations gathered in the following dataset: ((ad, y l ), i — 1,..., n). 
Ultimately Kriging aims at producing a predictor (also known as a response surface or 
metamodel) denoted as G. 

It is assumed that the model G is a realization of the normal process Y : i? x R d —» M 
defined by: 

Y(oj, x) — m(x) + Z(oj, x) (12) 

where m(x) is the trend and Z(x) is a zero-mean Gaussian process with a covariance 
function cq : x f d —y I which depends on the vector of parameters 6 G MY 6 : 

Z(y )] = cq(x, y ) (13) 

The trend is generally taken equal to the generalized linear model: 

m(x) = (f(x)) t p (14) 

where (/(a?)) f = and (3 = (fii,, /3 P ). Then, the Kriging method approxi¬ 

mates the model / by the mean of the Y given that: 

Y(oj,x^)—y^ Vi = l,...,n (15) 

The Kriging meta model G of G writes: 

G(x) = E[Y(u;, x)\Y'(u, x^) = y^ l \ Vi — 1,... ,n] (16) 

The meta model is then defined by: 

G(x) = (f{x)) t ^+(c e (x)) t C^ 1 {y-F^) (17) 

where 0 is the least squares estimator for f3 defined by: 
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E>={F , Cg 1 F) l F‘Cg 1 y (18) 

and Cq = [cg(xi,Xj)]i tj= i... n , F = [f{x i ) t } i=L .. n and c^(a;) = [cQ{x,Xi)\ i= The line 
command writes: 


»> algo = KrigingAlgorithm (inputSample , outputSample , basis, covarianceModel) 
»> algo, run () 

»> result = algo . getResult () 

»> metamodel = result . getMetaModel () 

»> graph = metamodel. draw () 


Figure 19 approximates the previously defined model Q : x i —> x sin x with a 


realization of a Gaussian process based on 6 observations. 



Fig. 19. An example of kriging approximation based on 6 observations: the blue line is the reference 
function Q : x <—> a; sin a; and the red one its approximation by a realization of a Gaussian process. 






Using an external code 
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The external simulator 

Fast evaluations of G 

On the practical side, the OpenTURNS software provides features which make the 
connection to the simulator G easy, and make its evaluation generally fast. Within the 
OpenTURNS framework, the method to connect to G is called "wrapping”. 

In the simplest situations, the function G is analytical and the formulas can be 
provided to OpenTURNS with a character string. Here, the Muparser C++ library [I] 
is used to evaluate the value of the mathematical function. In this case, the evaluation 
of G by OpenTURNS is quite fast. 

In the following Python script, consider the function G : M 3 —> M 2 , where 
G i(x) = x\ + X2 + X3 and G 2 (x) = x\ — X2X3, for any real numbers aq, X2 and X3. 
The input argument of the NumericalMathFunction class is a Python tuple, where 
the first item describes the three inputs variables, the second item describes the two 
output variables and the last item describes the two functions G\ and G 2 . 

»>G = NumericalMathFunction ( 

("xO" , " xl" ,"x2" ) , 

(” yO ” ,"yl" ) , 

("x0+xl+x2",”xO—xl*x2")) 

Once created, the function G can be used as a regular Python function, or can be passed as 
an input argument of other OpenTURNS classes. 

In most cases, the function G is provided as a Python function, which can be connected to 
OpenTURNS with the PythonFunction class. This task is easy (for those who are familiar with 
this language), and allows the scientific packages already available in Python to be combined. For 
example, if the computational code uses XML files on input or output, it is easy to make use of 
the XML features of Python (e.g. the minidom package). Moreover, if the function evaluation can be 
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vectorized (e.g. with the numpy package), then the func_sample option of the PythonFunction class 
can improve the performance a lot. 

The following Python script creates the function G associated with the flooding model. The 
flood function is first defined with the def Python statement. This function takes the variable X as 
input argument, which is an array with four components, Q, K_s, Z_v and Z_m, which corresponds to 
the input random variables in the model. The body of the flood function is a regular Python script, 
so that all Python functions can be used at this point (e.g. the numpy or scipy functions). The last 
statement of the function returns the overflow S. Then the PythonFunction class is used in order 
to convert this Python function into an object that OpenTURNS can use. This class takes as input 
arguments the number of input variables (in this case, 4), the number of outputs (in this case, 1) and 
the function and returns the object G. 

»>from openturns import PythonFunction 
»>def flood (X) : 

L = 5.0 e3; B = 300.0 
Q, K_s, Z_v, Z_m = X 
alpha = (Z__m — Z_v)/L 
H = (Q/(K_s*B* sqrt (alpha )))**( 0.6) 
return H 

»>G = PythonFunction (4 , 1, flood) 

If, as many of the computational codes commonly used, the data exchange is based on text 
files, OpenTURNS provides a component (coupling_tools) which is able to read and write structured 
text files based, for example, on line indices and perhaps containing tables (using line and column 
indices). Moreover, OpenTURNS provides a component which can evaluate such a Python function 
using the multi-thread capabilities that most computers have. 

Finally, when the computational code G is provided as a C or Fortran library, OpenTURNS 
provides a generic component to exchange data by memory, which is much faster than with files. In 
this case, the evaluation of G is automatically multi-thread. This component can be configured by 
Python, based on a XML file. If this method is not flexible enough, then the connection can be done 
with the C++ library. 

The previous techniques are documented in the OpenTURNS Developer’s Guide [T|. 
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Figure [20] compares the performance of three methods to connect to the function G : the 
PythonFunction class, the PythonFunction class with the func_sample option and the analytical 
function. This test was performed with a basic MS Windows laptop computer. Obviously, the fastest 
method is the analytical function, which can provide as many as 0.2 million samples per second, a 
performance which is four times the performance of the PythonFunction class. 



Fig. 20. Performance of various connection methods in OpenTURNS 


Evaluation of the derivatives 

When the algorithm involves an optimization step (e.g. in the FORM-SORM method) or a local 
approximation of G (e.g. in the Taylor development used to approximate the expectation and variance), 
the derivatives of G are required. 

When the computer code can compute the gradient and Hessian matrix of G, this information 
can be used by OpenTURNS . This happens sometimes, for example when the computer code has 
been differentiated with automatic differentiation methods, such as forward or adjoint techniques. 

In the case where the function is analytical and is provided as a character string, OpenTURNS 
is able to compute the exact derivatives of G. In order to do this, the software uses the Ev3 C++ 
library [2Sj to perform the symbolic computation of the derivatives and MuParser [2] to evaluate it. 
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In most common situations, however, the code does not compute its derivatives. In this case, 
OpenTURNS provides a method to compute the derivatives based on finite difference formulas. By 
default a centered finite difference formula for the gradient and a centered formula for the Hessian 
matrix are used. 

High performance computing 

For most common engineering practices, OpenTURNS can evaluate G with the multi-thread capa¬ 
bilities of most laptop and scientific workstations. However, when the evaluation of G is more CPU 
consuming or when the number of evaluations required is larger, these features are not sufficient by 
themselves and it is necessary to use a high performance computer such as the Zumbrota, Athos or 
Ivanhoe supercomputers available at EDF R&D which have from 16 000 to 65 000 cores m- 

In this case, two solutions are commonly used. The first one is to use a feature which can 
execute a Python function on remote processors, connected on the network with ssh. Here, the data 
flow is based on files, located in automatically generated directories, which prevents the loss of inter¬ 
mediate data. This feature (the DistributedPythonFunction) allows each remote processor to use 
its multi-thread capabilities, providing two different levels of parallel computing. 

The second solution is to use the OpenTURNS component integrated in the Salome platform. 
This component, along with a graphical user interface, called "Eficas", makes use of a software, called 
"YACS", which can call a Python script. The YACS module allows calculation schemes in Salome to 
be built, edited and executed. It provides both a graphical user interface to chain the computations by 
linking the inputs and outputs of computer codes, and then to execute these computations on remote 
machines. 

Several studies have been conducted at EDF based on the OpenTURNS component of Salome. 
For example, an uncertainty propagation study (the thermal evaluation of the storage of high-level 
nuclear waste) was making use of a computer code where one single run required approximately 10 
minutes on the 8 cores of a workstation (with shared memory). Within Salome, the OpenTURNS 
simulation involving 6000 unitary evaluations of the function G required 8000 CPU hours on 32 nodes 
0 - 


Conclusions 


42 


This educational example has shown a number of questions and problems that can be addressed by 
UQ methods: uncertainty quantification, central tendency evaluation, excess probability assessment 
and sensitivity analysis, that can require the use of a metamodel. 

Different numerical methods have been used for solving these three classes of problems, lead¬ 
ing substantially to the same (or very similar) results. In the industrial practice of UQ, the main issue 
(which actually motivates the choice of one mathematical method instead of another) is the compu¬ 
tational budget, which is actually given by the number of allowed runs of the deterministic model 
G. When the computer code implementing G is computationally expensive, one needs specifically 
designed mathematical and software tools. 

OpenTURNS is specially intended to meet these issues : (i) it includes a set of efficient 
mathematical methods for UQ and (ii) it can be easily connected to any external black box model 
G. Thanks to these two main features, OpenTURNS is a software that can address many different 
physics problems, and thus help to solve industrial problems. From this perspective, the partnership 
around OpenTURNS focuses efforts on the integration of the most efficient and innovative methods 
required by the industrial applications that takes into account both the need of genericity and of ease 
to communicate. The main projects for 2015 concern the improvement of the kriging implementa¬ 
tion to integrate some very smart methods of optimization. Around this theme some other classical 
optimization methods will also be generalized or newly implemented. 

A growing need for model exploration and analysis of uncertainty problem in industrial appli¬ 
cations is to better visualise the information provided by such a volume of data. In this area, specific 
visualization software, such as Paraview, can provide very efficient and interactive features. Taking 
the advantage of the integration of OpenTURNS in the Salome platform, EDF is working on a bet¬ 
ter link between the Paraview module in Salome (called Para VIS) and the uncertainty analysis with 
OpenTURNS : in 2012, functional boxplot m) has been implemented. Some recent work around 
in-situ visualization for uncertainty analysis should also be developed and implemented and so benefit 
very computationaly expensive model physics that generate an extremely high volume of data. 

Part of this work has been backed by French National Research Agency (ANR) trough the 
Chorus project (no. ANR-13-MONU-0005-08). We are grateful to the OpenTURNS Consortium mem¬ 
bers. We also thank Regis Lebrun, Mathieu Couplet and Merlin Keller for their help. 
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